#include <iostream>
#include <cmath>
#include <string>
#include <iomanip>
using namespace std;
double* devidedDifference(double* x, double* y,  int n) {
    int k = 0;
    double* f=new double [n*(n-1)/2];
    for (int i = 0; i < n - 1; i++) {
        for (int j = 0; j < n - i - 1; j++) {
            if (i == 0) {
                f[k] = (y[j + 1] - y[j]) / (x[j + 1] - x[j]);
            } else {
                f[k] = (f[k - (n - i - 1)] - f[k - (n - i - 1) - 1]) / (x[j + i + 1] - x[j]);
            }
            k++;
        }
    }
    return f;
}
double Newton(double* x, double* y, int n, double a) {
    double* b = new double[n];
    double* f=new double [n*(n-1)/2];
    f=devidedDifference(x,y,n);
    b[0] = 1;
    for (int i = 0; i < n - 1; i++) {
        b[i + 1] = b[i] * (a - x[i]);
    }
    double result = y[0];
    int k = 0;
    for (int i = 1; i < n; i++) {
        result += b[i] * f[k];
        k += (n - i);
    }
    delete[] f;
    delete[] b; // 释放动态分配的内存
    return result;
}

string newtonPolynomialExpression(double* x, double* y,int n) {
    double* f = devidedDifference(x, y, n);
    string expression = "";

    for (int i = 0; i < n; ++i) {
        if (i > 0) {
            expression += " + ";
        }
        expression += to_string(f[i]);

        for (int j = 0; j < i; ++j) {
            expression += "(a - " + to_string(x[j]) + ")";
        }
    }

    return expression;
}
int main(){
    int n=7;
    double* x=new double[n];
    double* y1=new double[n];
    double* y2=new double[n];
    double init_x[] = {0, 6,10,13,17,20,28};
    copy(init_x, init_x + n, x);

    double init_y1[] = {6.67,17.3,42.7,37.3,30.1,29.3,28.7};
    copy(init_y1, init_y1 + n, y1);

    double init_y2[] = {6.67,16.7,18.9,15.0,10.6,9.44,8.89};
    copy(init_y2, init_y2 + n, y2);
    cout<<"(a)"<<endl;
    setprecision(2);
    cout<<"The polynomial for Sp1: "<<newtonPolynomialExpression(x,y1,n)<<endl;
    cout<<"The polynomial for Sp2: "<<newtonPolynomialExpression(x,y2,n)<<endl;
    cout<<"(b)"<<endl;
    //a=43
    int a=43;
    double result1=Newton(x,y1,n,a);
    double result2=Newton(x,y2,n,a);
    if(result1>0){
        if(result2>0){
            cout<<"Neither sample of larvae will die after another 15 days."<<endl;
        }else{
            cout<<"The second sample will die after another 15 days."<<endl;
        }
    }else{
        cout<<"Both sample of larvae will die after another 15 days."<<endl;
    }
    return 0;
}
